//---------------------------------------------------------------------------------------------------
template <class Real>
Matrix4<Real>::Matrix4(bool bZero)
{
	if (bZero)
	{
		MakeZero();
	} 
	else
	{
		MakeIdentity();
	}
}
//---------------------------------------------------------------------------------------------------
template <class Real>
Matrix4<Real>::Matrix4(Real f00, Real f01, Real f02, Real f03,
				 Real f10, Real f11, Real f12, Real f13,
				 Real f20, Real f21, Real f22, Real f23,
				 Real f30, Real f31, Real f32, Real f33)
{
	m_afTuple[0] = f00;
	m_afTuple[1] = f01;
	m_afTuple[2] = f02;
	m_afTuple[3] = f03;
	m_afTuple[4] = f10;
	m_afTuple[5] = f11;
	m_afTuple[6] = f12;
	m_afTuple[7] = f13;
	m_afTuple[8] = f20;
	m_afTuple[9] = f21;
	m_afTuple[10] = f22;
	m_afTuple[11] = f23;
	m_afTuple[12] = f30;
	m_afTuple[13] = f31;
	m_afTuple[14] = f32;
	m_afTuple[15] = f33;
}
//---------------------------------------------------------------------------------------------------
template <class Real>
Matrix4<Real>::Matrix4(const Matrix4& rkM)
{
	m_afTuple[0] = rkM.m_afTuple[0];
	m_afTuple[1] = rkM.m_afTuple[1];
	m_afTuple[2] = rkM.m_afTuple[2];
	m_afTuple[3] = rkM.m_afTuple[3];
	m_afTuple[4] = rkM.m_afTuple[4];
	m_afTuple[5] = rkM.m_afTuple[5];
	m_afTuple[6] = rkM.m_afTuple[6];
	m_afTuple[7] = rkM.m_afTuple[7];
	m_afTuple[8] = rkM.m_afTuple[8];
	m_afTuple[9] = rkM.m_afTuple[9];
	m_afTuple[10] = rkM.m_afTuple[10];
	m_afTuple[11] = rkM.m_afTuple[11];
	m_afTuple[12] = rkM.m_afTuple[12];
	m_afTuple[13] = rkM.m_afTuple[13];
	m_afTuple[14] = rkM.m_afTuple[14];	
	m_afTuple[15] = rkM.m_afTuple[15];
}
//---------------------------------------------------------------------------------------------------
template <class Real>
Matrix4<Real>::Matrix4(const Vector3<Real>& rkAxis, Real fAngleRad)
{
	FromAxisAngle(rkAxis, fAngleRad);
}
//---------------------------------------------------------------------------------------------------
template <class Real>
Matrix4<Real>& Matrix4<Real>::MakeZero()
{
	m_afTuple[0] = 0.0f;
	m_afTuple[1] = 0.0f;
	m_afTuple[2] = 0.0f;
	m_afTuple[3] = 0.0f;
	m_afTuple[4] = 0.0f;
	m_afTuple[5] = 0.0f;
	m_afTuple[6] = 0.0f;
	m_afTuple[7] = 0.0f;
	m_afTuple[8] = 0.0f;
	m_afTuple[9] = 0.0f;
	m_afTuple[10] = 0.0f;
	m_afTuple[11] = 0.0f;
	m_afTuple[12] = 0.0f;
	m_afTuple[13] = 0.0f;
	m_afTuple[14] = 0.0f;
	m_afTuple[15] = 0.0f;
	return *this;
}
//---------------------------------------------------------------------------------------------------
template <class Real>
Matrix4<Real>& Matrix4<Real>::MakeIdentity()
{
	m_afTuple[0] = 1.0f;
	m_afTuple[1] = 0.0f;
	m_afTuple[2] = 0.0f;
	m_afTuple[3] = 0.0f;
	m_afTuple[4] = 0.0f;
	m_afTuple[5] = 1.0f;
	m_afTuple[6] = 0.0f;
	m_afTuple[7] = 0.0f;
	m_afTuple[8] = 0.0f;
	m_afTuple[9] = 0.0f;
	m_afTuple[10] = 1.0f;
	m_afTuple[11] = 0.0f;
	m_afTuple[12] = 0.0f;
	m_afTuple[13] = 0.0f;
	m_afTuple[14] = 0.0f;
	m_afTuple[15] = 1.0f;
	return *this;
}
//---------------------------------------------------------------------------------------------------
template <class Real>
int Matrix4<Real>::CompareArrays(const Matrix4& rkM) const
{
	return memcmp(m_afTuple, rkM.m_afTuple, 16*sizeof(Real));
}
//---------------------------------------------------------------------------------------------------
template <class Real>
bool Matrix4<Real>::operator== (const Matrix4& rkM) const
{
	return CompareArrays(rkM) == 0;
}
//---------------------------------------------------------------------------------------------------
template <class Real>
bool Matrix4<Real>::operator!= (const Matrix4& rkM) const
{
	return CompareArrays(rkM) != 0;
}
//---------------------------------------------------------------------------------------------------
template <class Real>
bool Matrix4<Real>::operator< (const Matrix4& rkM) const
{
	return CompareArrays(rkM) < 0;
}
//---------------------------------------------------------------------------------------------------
template <class Real>
bool Matrix4<Real>::operator<= (const Matrix4& rkM) const
{
	return CompareArrays(rkM) <= 0;
}
//---------------------------------------------------------------------------------------------------
template <class Real>
bool Matrix4<Real>::operator> (const Matrix4& rkM) const
{
	return CompareArrays(rkM) > 0;
}
//---------------------------------------------------------------------------------------------------
template <class Real>
bool Matrix4<Real>::operator>= (const Matrix4& rkM) const
{
	return CompareArrays(rkM) >= 0;
}
//---------------------------------------------------------------------------------------------------
template <class Real>
Matrix4<Real>& Matrix4<Real>::FromAxisAngle(const Vector3<Real>& rkAxis, Real fAngleRad)
{
	Real fCos = (Real)cos((double)fAngleRad);
	Real fSin = (Real)sin((double)fAngleRad);
	Real fOneMinusCos = 1.0f - fCos;
	Real fX2 = rkAxis[0]*rkAxis[0];
	Real fY2 = rkAxis[1]*rkAxis[1];
	Real fZ2 = rkAxis[2]*rkAxis[2];
	Real fXYM = rkAxis[0]*rkAxis[1]*fOneMinusCos;
	Real fXZM = rkAxis[0]*rkAxis[2]*fOneMinusCos;
	Real fYZM = rkAxis[1]*rkAxis[2]*fOneMinusCos;
	Real fXSin = rkAxis[0]*fSin;
	Real fYSin = rkAxis[1]*fSin;
	Real fZSin = rkAxis[2]*fSin;	

	m_afTuple[ 0] = fX2*fOneMinusCos + fCos;
	m_afTuple[ 1] = fXYM + fZSin;
	m_afTuple[ 2] = fXZM - fYSin;
	m_afTuple[ 3] = 0.0f;

	m_afTuple[ 4] = fXYM - fZSin;
	m_afTuple[ 5] = fY2*fOneMinusCos + fCos;
	m_afTuple[ 6] = fYZM + fXSin;
	m_afTuple[ 7] = 0.0f;

	m_afTuple[ 8] = fXZM + fYSin;
	m_afTuple[ 9] = fYZM - fXSin;
	m_afTuple[10] = fZ2*fOneMinusCos + fCos;
	m_afTuple[11] = 0.0f;

	m_afTuple[12] = 0.0f;
	m_afTuple[13] = 0.0f;
	m_afTuple[14] = 0.0f;
	m_afTuple[15] = 1.0f;

	return *this;
}
//---------------------------------------------------------------------------------------------------
template <class Real>
void Matrix4<Real>::Print()
{
	for (int i = 1; i<= 16; i++)
	{
		printf("%6f ", m_afTuple[i-1]);
		if (i%4 == 0)
		{
			printf("\n");
		}
	}
	printf("\n");
}
//---------------------------------------------------------------------------------------------------
template <class Real>
Matrix4<Real>& Matrix4<Real>::MakeHomogeneous(Matrix3<Real>& rkM3)
{
	m_afTuple[0] = rkM3[0];
	m_afTuple[1] = rkM3[1];
	m_afTuple[2] = rkM3[2];
	m_afTuple[3] = 0.0f;

	m_afTuple[4] = rkM3[3];
	m_afTuple[5] = rkM3[4];
	m_afTuple[6] = rkM3[5];
	m_afTuple[7] = 0.0f;

	m_afTuple[8] = rkM3[6];
	m_afTuple[9] = rkM3[7];
	m_afTuple[10] = rkM3[8];
	m_afTuple[11] = 0.0f;

	m_afTuple[12] = 0.0f;
	m_afTuple[13] = 0.0f;
	m_afTuple[14] = 0.0f;
	m_afTuple[15] = 1.0f;

	return *this;
}
//---------------------------------------------------------------------------------------------------
//---------------------------------------------------------------------------------------------------
template <class Real>
inline Real& Matrix4<Real>::operator [] (int iEntry)
{
	if (iEntry < 0)
	{
		iEntry = 0;
	}
	if (iEntry > 15)
	{
		iEntry = 15;
	}
	return m_afTuple[iEntry];
}
//---------------------------------------------------------------------------------------------------
template <class Real>
inline Matrix4<Real>::operator const Real* () const
{
	return m_afTuple;
}
//---------------------------------------------------------------------------------------------------
template <class Real>
inline Matrix4<Real>::operator Real* ()
{
	return m_afTuple;
}
//---------------------------------------------------------------------------------------------------
template <class Real>
inline Matrix4<Real> Matrix4<Real>::operator* (const Matrix4& rkM) const
{
	return Matrix4(
		m_afTuple[ 0]*rkM.m_afTuple[ 0] +
		m_afTuple[ 1]*rkM.m_afTuple[ 4] +
		m_afTuple[ 2]*rkM.m_afTuple[ 8] +
		m_afTuple[ 3]*rkM.m_afTuple[12],

		m_afTuple[ 0]*rkM.m_afTuple[ 1] +
		m_afTuple[ 1]*rkM.m_afTuple[ 5] +
		m_afTuple[ 2]*rkM.m_afTuple[ 9] +
		m_afTuple[ 3]*rkM.m_afTuple[13],

		m_afTuple[ 0]*rkM.m_afTuple[ 2] +
		m_afTuple[ 1]*rkM.m_afTuple[ 6] +
		m_afTuple[ 2]*rkM.m_afTuple[10] +
		m_afTuple[ 3]*rkM.m_afTuple[14],

		m_afTuple[ 0]*rkM.m_afTuple[ 3] +
		m_afTuple[ 1]*rkM.m_afTuple[ 7] +
		m_afTuple[ 2]*rkM.m_afTuple[11] +
		m_afTuple[ 3]*rkM.m_afTuple[15],

		m_afTuple[ 4]*rkM.m_afTuple[ 0] +
		m_afTuple[ 5]*rkM.m_afTuple[ 4] +
		m_afTuple[ 6]*rkM.m_afTuple[ 8] +
		m_afTuple[ 7]*rkM.m_afTuple[12],

		m_afTuple[ 4]*rkM.m_afTuple[ 1] +
		m_afTuple[ 5]*rkM.m_afTuple[ 5] +
		m_afTuple[ 6]*rkM.m_afTuple[ 9] +
		m_afTuple[ 7]*rkM.m_afTuple[13],

		m_afTuple[ 4]*rkM.m_afTuple[ 2] +
		m_afTuple[ 5]*rkM.m_afTuple[ 6] +
		m_afTuple[ 6]*rkM.m_afTuple[10] +
		m_afTuple[ 7]*rkM.m_afTuple[14],

		m_afTuple[ 4]*rkM.m_afTuple[ 3] +
		m_afTuple[ 5]*rkM.m_afTuple[ 7] +
		m_afTuple[ 6]*rkM.m_afTuple[11] +
		m_afTuple[ 7]*rkM.m_afTuple[15],

		m_afTuple[ 8]*rkM.m_afTuple[ 0] +
		m_afTuple[ 9]*rkM.m_afTuple[ 4] +
		m_afTuple[10]*rkM.m_afTuple[ 8] +
		m_afTuple[11]*rkM.m_afTuple[12],

		m_afTuple[ 8]*rkM.m_afTuple[ 1] +
		m_afTuple[ 9]*rkM.m_afTuple[ 5] +
		m_afTuple[10]*rkM.m_afTuple[ 9] +
		m_afTuple[11]*rkM.m_afTuple[13],

		m_afTuple[ 8]*rkM.m_afTuple[ 2] +
		m_afTuple[ 9]*rkM.m_afTuple[ 6] +
		m_afTuple[10]*rkM.m_afTuple[10] +
		m_afTuple[11]*rkM.m_afTuple[14],

		m_afTuple[ 8]*rkM.m_afTuple[ 3] +
		m_afTuple[ 9]*rkM.m_afTuple[ 7] +
		m_afTuple[10]*rkM.m_afTuple[11] +
		m_afTuple[11]*rkM.m_afTuple[15],

		m_afTuple[12]*rkM.m_afTuple[ 0] +
		m_afTuple[13]*rkM.m_afTuple[ 4] +
		m_afTuple[14]*rkM.m_afTuple[ 8] +
		m_afTuple[15]*rkM.m_afTuple[12],

		m_afTuple[12]*rkM.m_afTuple[ 1] +
		m_afTuple[13]*rkM.m_afTuple[ 5] +
		m_afTuple[14]*rkM.m_afTuple[ 9] +
		m_afTuple[15]*rkM.m_afTuple[13],

		m_afTuple[12]*rkM.m_afTuple[ 2] +
		m_afTuple[13]*rkM.m_afTuple[ 6] +
		m_afTuple[14]*rkM.m_afTuple[10] +
		m_afTuple[15]*rkM.m_afTuple[14],

		m_afTuple[12]*rkM.m_afTuple[ 3] +
		m_afTuple[13]*rkM.m_afTuple[ 7] +
		m_afTuple[14]*rkM.m_afTuple[11] +
		m_afTuple[15]*rkM.m_afTuple[15]);
}
//---------------------------------------------------------------------------------------------------
template <class Real>
inline Matrix4<Real> Matrix4<Real>::operator* (Real fScalar) const
{
	return Matrix4(
		fScalar*m_afTuple[0],
		fScalar*m_afTuple[1],
		fScalar*m_afTuple[2],
		fScalar*m_afTuple[3],
		fScalar*m_afTuple[4],
		fScalar*m_afTuple[5],
		fScalar*m_afTuple[6],
		fScalar*m_afTuple[7],
		fScalar*m_afTuple[8],
		fScalar*m_afTuple[9],
		fScalar*m_afTuple[10],
		fScalar*m_afTuple[11],
		fScalar*m_afTuple[12],
		fScalar*m_afTuple[13],
		fScalar*m_afTuple[14],
		fScalar*m_afTuple[15]);
}
//---------------------------------------------------------------------------------------------------
template <class Real>
inline Vector4<Real> Matrix4<Real>::operator* (const Vector4<Real>& rkV) const
{
	return Vector4<Real>(
		m_afTuple[ 0]*rkV[0] +
		m_afTuple[ 4]*rkV[1] +
		m_afTuple[ 8]*rkV[2] +
		m_afTuple[12]*rkV[3],

		m_afTuple[ 1]*rkV[0] +
		m_afTuple[ 5]*rkV[1] +
		m_afTuple[ 9]*rkV[2] +
		m_afTuple[13]*rkV[3],

		m_afTuple[ 2]*rkV[0] +
		m_afTuple[ 6]*rkV[1] +
		m_afTuple[10]*rkV[2] +
		m_afTuple[14]*rkV[3],

		m_afTuple[ 3]*rkV[0] +
		m_afTuple[ 7]*rkV[1] +
		m_afTuple[11]*rkV[2] +
		m_afTuple[15]*rkV[3]);
}
//---------------------------------------------------------------------------------------------------
template <class Real>
inline Matrix4<Real> Matrix4<Real>::Transpose() const
{
	return Matrix4(
		m_afTuple[ 0],
		m_afTuple[ 4],
		m_afTuple[ 8],
		m_afTuple[12],
		m_afTuple[ 1],
		m_afTuple[ 5],
		m_afTuple[ 9],
		m_afTuple[13],
		m_afTuple[ 2],
		m_afTuple[ 6],
		m_afTuple[10],
		m_afTuple[14],
		m_afTuple[ 3],
		m_afTuple[ 7],
		m_afTuple[11],
		m_afTuple[15]);
}
//---------------------------------------------------------------------------------------------------
template <class Real>
inline Matrix4<Real> Matrix4<Real>::Inverse() const
{
    Real fA0 = m_afTuple[ 0]*m_afTuple[ 5] - m_afTuple[ 1]*m_afTuple[ 4];
    Real fA1 = m_afTuple[ 0]*m_afTuple[ 6] - m_afTuple[ 2]*m_afTuple[ 4];
    Real fA2 = m_afTuple[ 0]*m_afTuple[ 7] - m_afTuple[ 3]*m_afTuple[ 4];
    Real fA3 = m_afTuple[ 1]*m_afTuple[ 6] - m_afTuple[ 2]*m_afTuple[ 5];
    Real fA4 = m_afTuple[ 1]*m_afTuple[ 7] - m_afTuple[ 3]*m_afTuple[ 5];
    Real fA5 = m_afTuple[ 2]*m_afTuple[ 7] - m_afTuple[ 3]*m_afTuple[ 6];
    Real fB0 = m_afTuple[ 8]*m_afTuple[13] - m_afTuple[ 9]*m_afTuple[12];
    Real fB1 = m_afTuple[ 8]*m_afTuple[14] - m_afTuple[10]*m_afTuple[12];
    Real fB2 = m_afTuple[ 8]*m_afTuple[15] - m_afTuple[11]*m_afTuple[12];
    Real fB3 = m_afTuple[ 9]*m_afTuple[14] - m_afTuple[10]*m_afTuple[13];
    Real fB4 = m_afTuple[ 9]*m_afTuple[15] - m_afTuple[11]*m_afTuple[13];
    Real fB5 = m_afTuple[10]*m_afTuple[15] - m_afTuple[11]*m_afTuple[14];

    Real fDet = fA0*fB5-fA1*fB4+fA2*fB3+fA3*fB2-fA4*fB1+fA5*fB0;
    if ((Real)fabs(fDet) <= Math<Real>::ZERO_TOLERANCE)
    {
        return Matrix4::ZERO;
    }

    Matrix4 kInv;
    kInv.m_afTuple[ 0] =
        + m_afTuple[ 5]*fB5 - m_afTuple[ 6]*fB4 + m_afTuple[ 7]*fB3;
    kInv.m_afTuple[ 4] =
        - m_afTuple[ 4]*fB5 + m_afTuple[ 6]*fB2 - m_afTuple[ 7]*fB1;
    kInv.m_afTuple[ 8] =
        + m_afTuple[ 4]*fB4 - m_afTuple[ 5]*fB2 + m_afTuple[ 7]*fB0;
    kInv.m_afTuple[12] =
        - m_afTuple[ 4]*fB3 + m_afTuple[ 5]*fB1 - m_afTuple[ 6]*fB0;
    kInv.m_afTuple[ 1] =
        - m_afTuple[ 1]*fB5 + m_afTuple[ 2]*fB4 - m_afTuple[ 3]*fB3;
    kInv.m_afTuple[ 5] =
        + m_afTuple[ 0]*fB5 - m_afTuple[ 2]*fB2 + m_afTuple[ 3]*fB1;
    kInv.m_afTuple[ 9] =
        - m_afTuple[ 0]*fB4 + m_afTuple[ 1]*fB2 - m_afTuple[ 3]*fB0;
    kInv.m_afTuple[13] =
        + m_afTuple[ 0]*fB3 - m_afTuple[ 1]*fB1 + m_afTuple[ 2]*fB0;
    kInv.m_afTuple[ 2] =
        + m_afTuple[13]*fA5 - m_afTuple[14]*fA4 + m_afTuple[15]*fA3;
    kInv.m_afTuple[ 6] =
        - m_afTuple[12]*fA5 + m_afTuple[14]*fA2 - m_afTuple[15]*fA1;
    kInv.m_afTuple[10] =
        + m_afTuple[12]*fA4 - m_afTuple[13]*fA2 + m_afTuple[15]*fA0;
    kInv.m_afTuple[14] =
        - m_afTuple[12]*fA3 + m_afTuple[13]*fA1 - m_afTuple[14]*fA0;
    kInv.m_afTuple[ 3] =
        - m_afTuple[ 9]*fA5 + m_afTuple[10]*fA4 - m_afTuple[11]*fA3;
    kInv.m_afTuple[ 7] =
        + m_afTuple[ 8]*fA5 - m_afTuple[10]*fA2 + m_afTuple[11]*fA1;
    kInv.m_afTuple[11] =
        - m_afTuple[ 8]*fA4 + m_afTuple[ 9]*fA2 - m_afTuple[11]*fA0;
    kInv.m_afTuple[15] =
        + m_afTuple[ 8]*fA3 - m_afTuple[ 9]*fA1 + m_afTuple[10]*fA0;

    Real fInvDet = ((Real)1.0)/fDet;
    kInv.m_afTuple[ 0] *= fInvDet;
    kInv.m_afTuple[ 1] *= fInvDet;
    kInv.m_afTuple[ 2] *= fInvDet;
    kInv.m_afTuple[ 3] *= fInvDet;
    kInv.m_afTuple[ 4] *= fInvDet;
    kInv.m_afTuple[ 5] *= fInvDet;
    kInv.m_afTuple[ 6] *= fInvDet;
    kInv.m_afTuple[ 7] *= fInvDet;
    kInv.m_afTuple[ 8] *= fInvDet;
    kInv.m_afTuple[ 9] *= fInvDet;
    kInv.m_afTuple[10] *= fInvDet;
    kInv.m_afTuple[11] *= fInvDet;
    kInv.m_afTuple[12] *= fInvDet;
    kInv.m_afTuple[13] *= fInvDet;
    kInv.m_afTuple[14] *= fInvDet;
    kInv.m_afTuple[15] *= fInvDet;

    return kInv;
}
//---------------------------------------------------------------------------------------------------
template <class Real>
inline Vector4<Real> operator* (const Vector4<Real>& rkV,
    const Matrix4<Real>& rkM)
{
	return Vector4<Real>(
		rkV[0]*rkM[0]+rkV[1]*rkM[4]+rkV[2]*rkM[8]+rkV[3]*rkM[12],
		rkV[0]*rkM[1]+rkV[1]*rkM[5]+rkV[2]*rkM[9]+rkV[3]*rkM[13],
		rkV[0]*rkM[2]+rkV[1]*rkM[6]+rkV[2]*rkM[10]+rkV[3]*rkM[14],
		rkV[0]*rkM[3]+rkV[1]*rkM[7]+rkV[2]*rkM[11]+rkV[3]*rkM[15]);
}
//---------------------------------------------------------------------------------------------------